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Oscillation is an important cellular process that regulates timing of different vital 
life cycles. However, in the noisy cellular environment, oscillations can be highly 
inaccurate due to phase fluctuations. It remains poorly understood how biochemical 
circuits suppress phase fluctuations and what is the incurred thermodynamic cost. 
Here, we study four different types of biochemical oscillations representing three 
basic oscillation motifs shared by all known oscillatory systems. We find that the 
phase diffusion constant follows the same inverse dependence on the free energy 
dissipation per period for all systems studied. This relationship between the phase 
diffusion and energy dissipation is shown analytically in a model of noisy oscillation. 
Microscopically, we hnd that the oscillation is driven by multiple irreversible cycles 
that hydrolyze the fuel molecules such as ATP; the number of phase coherent periods 
is proportional to the free energy consumed per period. Experimental evidence in 
support of this universal relationship and testable predictions are also presented. 
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I. INTRODUCTION 


Living systems are dissipative, consuming ener^ to perform key functions for their sur¬ 
vival and growth. While it is clear that free energy JlJ-yl is needed for physical functions, such 
as cell motility jd] and macromolecule synthesis js], it remains poorly understood whether 
and how regulatory functions are enhanced by free energy consumption. The relationship 
between biological regulatory functions and nonequilibrium thermodynamics has been an 
active area in biophysics jb-H]. For example, recent studies in different cellular adaptation 
processes demonstrated that the cost-performance trade-off follows a universal relationship, 
independent of the detailed biochemical circuits BB. 


Oscillatory behaviors exist in many biological systems, e.g., glycolysis 12|], cyclic AM 


signaling 13|, cell cycle 14j-ll6|. circadian rhythms (l2|,ll7|, and synthetic oscillators 18l Il9|. 


These biochemical oscillations are crucial in controlling the timing of life processes. Much 
is known now about the structure of biochemical circuits responsible for these oscillatory 
behaviors. There are a few basic network motifs, illustrated in Figure la, which are respon¬ 


sible for all known biochemical and genetic oscillations 


[ 3 ,[ 3 , 


16 


18| . These network motifs 


share a few essential features, such as nonlinearity, negative feedback, and a time delay, as 


summarized by Novak and Tyson in 


20| . However, in small systems such as a single cell, the 


dynamics of oscillations are subject to large fluctuations from the environment, due to their 
small sizes. Thus, one may ask how biological systems maintain coherence of oscillations 
amidst these fluctuations |2l|. Here, we study the thermodynamic cost of controlling oscil¬ 
lation coherence in different representative oscillatory systems and investigate whether there 
is a general (universal) relation between the accuracy of the oscillation and its minimum 
free energy cost that may apply to all biochemical oscillations. 


We study four specific models, the activator-inhibitor (AI) model, the repressilator model, 
the brusselator model, and the allosteric glycolysis model, chosen to exemplify the three 
different basic oscillation motifs, as shown in Fig. 1. For all the systems studied, a finite 
(critical) amount of free energy is needed to drive them to oscillate. Beyond the onset of 
oscillation, extra free energy dissipation is used to reduce the phase diffusion constant and 
thus enhance the coherence time and phase accuracy of the oscillations. A general inverse 
relationship between the phase diffusion constant and the free energy dissipation is found in 
all the four models studied, suggesting that the relation may hold true for all biochemical 
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oscillations. The energy-accuracy relation for noisy oscillations is also verified analytically 
in the noisy complex Stuart-Landau equation. In the following, we report these results 
followed by a in-depth discussion of a plausible general microscopic mechanism/strategy for 
energy-assisted noise suppression. 


II. MODELS AND RESULTS 


A. Four biochemical oscillators representing the basic network motifs 


All known biochemical and genetic oscillators contain at least one of the basic motifs (or 

n 

their variance) in network topology [201 ? ]. To search for general principles in these noisy 
oscillatory systems, we study four biochemical systems (Fig. [T]), each representing one of 
the three basic network motifs responsible for oscillatory behaviors. The first one is the 
activator-inhibitor (AI) system, where a negative feedback is interlinked with a positive 



this motif in a simplified biological network with a phosphorylation-dephosphorylation cycle 
(Fig. lb). The second model is a repressilator, which consists of three components connected 
in a negative feedback loop, such that each component represses the next one in the loop, 
and is itself repressed by the previous o ne ( Middle panel. Fig. la). The first synthetic ge¬ 


netic oscillator was built with this motif ISj]. Many important transcriptional-translational 
oscillators also use this motif as their backbone, such as circadian clock in mammalian cells 
[u|, NF-kB signaling 281, and the p53-mdm2 oscillations in cancer cells 291. Here, we take 
the repressilator composed of CDKl, Plkl, and APC in a cell cycle as our case study (Fig. 
Ic). The third model we chose is the brusselator, which is one of the simplest two-component 
systems that can generate sustained oscillations (Right panel, Fig, la). The concentration 
fluctuations due to small molecule numbers were analyzed in [SOj, here, we aim to study the 
effect of noise on the phase of the oscillation. The brusselator (Fig. Id) is a special kind of 
substrate-depletion model (3l|, where substrate S is converted by a process that is amplified 
autocatalytically by the product P. Exarnples of substrate-depletion motif are oscillations in 
glycolysis [l2, 32| and Calcium signaling j^. Here, we examine the noise effect in glycolysis 
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oscillations, where the allosteric enzyme PFK catalyzes snbstrate to prodnct in a network 
shown in Fig. le. 

In our study, we introduced a parameter 7 to characterize the reversibility of the biochem¬ 
ical networks. In a reaction loop, 7 corresponds to the ratio of the product of the reaction 
rates in one direction (e.g., counter-clock-wise) and that in the other direction (e.g., clock¬ 
wise). When 7 = 1, the system is in equilibrium without any free energy dissipation. For 
7 7 ^ 1, free energy is dissipated. Here, we study the relationship between the dynamics 
and the energetics of the biochemical networks by varying 7 . The mathematical details of 
the four models are described in the Supplemental Information (SI), all parameters (e.g., 
reaction rates, concentrations, time, and volumes) are shown here as dimensionless numbers 
with their units explained in SI. 


B. Phase diffusion reduces the coherence time 


In all four models that we studied, there is an onset of oscillation as 7 decreases below a 
critical value 7 c (< 1). This means that a hnite critical free energy dissipation [Wc > 0) is 
needed to generate an oscillatory behavior (see Fig. SI in SI). In Fig. 2a, two trajectories 
of the concentration of the inhibitor X are shown for 7 < 7 c in the activator-inhibitor 
model, where 7 c = 2 x 10“^. As evident in Fig. 2a, biochemical oscillations are noisy. 
To characterize the coherence of the oscillation in time, we computed the auto-correlation 
function C{t) for a given concentration variable x in the network. As shown in Fig. 2b, C{t) 
follows a damped oscillation: 


Cit) 


{{x{t + s) - (x))(x(g) - (x))), 
— (x)2 


exp(— t/uc) X cos(27rt/T), 


( 1 ) 


where T is the period and Tc dehnes a coherence time for the oscillation. 

The oscillatory state breaks time translation invariance (symmetry) of the underlying 
biochemical system. As a result, the phase of the oscillation is a soft mode and follows 
diffusive dynamics in the presence of noise. To quantify the phase diffusion, we simulated 
many trajectories in the model(s) with the same parameters and the same initial conditions. 
In Fig. 2c, the peak times for 500 trajectories in the AI model are shown in a raster plot 
together with the peak time distributions (red lines). The variance (a^) of the distribution 
versus the average peak time is shown in Fig. 2d. It is clear that the variance goes lin- 
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early with time, confirming the diffusive nature of the phase, and the linear slope dehnes a 
peak time diffusion constant D. It is easy to show that the coherence time Tc is inversely 
proportional to D: 

r, = «TVD, (2) 

where a is a constant dependent on the waveform {a = (27r^)“^ for a sine wave). 


C. Free energy dissipation suppresses phase diffusion 


As 7 decreases below 7 c, more free energy is dissipated. What is the effect of the additional 


free energy dissipation beyond the onset of oscil 
we can compute the free energy dissipation rate 


ation? From the chemical reaction rates. 


3 ^: 


w' = - ■’i 


In 


Jr 


( 3 ) 


where and are the forward and backward fluxes of the Ah reaction, and free energy 
is in units of ksT, set to unity here. For the activator-inhibitor and glycolysis models, we 
calculated the energy dissipation rate using Eq. [3l For systems with continuum stochas¬ 
tic dynamics described by Langevin equations (e.g., the brusselator and the repressilator 
models), we can obtain the steady-state distribution P{x) by solving the corresponding 
Fokker-Planck equation or by direct stochastic simulations (see Fig. S2 in SI for an ex¬ 
ample). From P{x), we computed the phase space fluxes and the free energy dissipation 
rate following {d] (see the Methods section and SI for details). For oscillatory systems, the 
dissipation rate W varies in a period T. We dehne AW = J^Wdt to characterize the free 
energy dissipation per period per volume. 

For each of the four models, AW and the dimensionless peak time diffusion constant D/T 
were computed for different parameter values (reaction rates, protein concentrations) in the 
oscillatory regime 7 < 7c and for different volume V. As shown in Fig. 3, for all the four 
models considered, D/T decreases as the energy dissipation AW increases and eventually 
saturates to a hxed value when AW —)■ 00 (i.e., 7 = 0 ). The phase diffusion constants scale 
inversely with the volume V. As shown in the insets of Fig. 3, the scaled D/T (by the 
volume V) collapsed onto a simple curve, which can be approximated by the same simple 
form in all the four models studied: 


V 

T 


AW-Wc’ 


( 4 ) 
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where Wc is the critical free energy, and Wq and C are intensive constants (independent of 
volume), whose values in different systems (models) are given in the legend of Fig. 3. 


D. The free energy sources and experimental evidence 


What is the free energy source driving the biochemical oscillations? For the activator- 
inhibitor model, the free energy is provided by ATP hydrolysis in the phosphorylation- 
dephosphorylation (PdP) cycle (see Fig. la). Besides the standard free energy AGq of 
ATP hydrolysis, the total free energy dissipation per period AhF also depends on (and 
thus can be controlled by) the concentrations of ATP, ADP and the inorganic phosphate 
Pj. These concentrations ([ATP], [ADP], and [Pj]) directly affect the biochemical reaction 
rates in our model and consequently the phase diffusion of the oscillation. In Fig. 4a, we 
show the phase diffusion constant {D/T) versus the dissipation per period (AW) for 300 
randomly chosen points in the oscillatory regime of the ([ATP], [ADP], [Pj]) space (see 
Fig. 4b). Remarkably, all the points lie above an envelope curve (the dotted line), which 
follows Eq|3l This envelope curve defines the best performance of the biochemical network, 
i.e., the minimum free energy AWm needed to achieve a given level of phase coherence. 
For each choice of the concentrations ([ATP], [APP], [PJ), a functional efficiency E can be 
dehned as the ratio of AWm and the actual cost AW for the same performance {D/T). The 
efficiency is represented by color in Fig. 4a&b. We investigated how efficiency depends on 
the three concentrations. As shown in Fig. 4c, the efficiency E does not simply increase 
with the ATP concentration; instead it peaks near a particular level of [ATP], at which the 
phosphorylation and dephosphorylation fluxes are matched. Similarly, E does not have any 
clear dependence on [APP] or [Pj] level, it is high near a fixed ratio of [APP]/[Pj], when 
the kinetic rates of the phosphorylation and dephosphorylation parts of the PdP cycle are 
matched. 

These predicted dependence of oscillatory behaviors on [ATP], [ADP], and [Pj] concen¬ 
trations, as shown in Fig. 4, may be tested experimentally by measuring peak-to-peak time 
variations or equivalently the correlation time for different nucleotide concentrations. As re¬ 
ported in two recent studies [35|, l36| , the oscillatory dynamics of the phosophorylated KaiC 
protein in a reconstituted circadian clock from cyanobacteria (the Kai system) have been 
measured in media with different ATP/ADP ratios. We analysed the data according to Eq. 
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1 and obtained the correlation time (tc) and the period (T) for different ATP/ADT ratios 
(see SI and Fig. S3 in SI for details). In Fig. [5l we plotted the period and the phase diffusion 
versus \n{[ATP]/[ADP]), which is the entropic contribution to the free energy dissipation. 
As the ATP/ADP ratio increases, the period changes little. In contrast, the phase diffusion 
T/tc = a~^D/T decreases significantly and eventually saturates at high ATP/ADP ratios, 
consistent with the relationship between energy dissipation and phase diffusion discovered 
here. 


E. Analytical results from the noisy Stuart-Landau equation 


To understand the relationship between phase accuracy and energy dissipation, we con¬ 
sider the noisy Stuart-Landau equation for a complex order parameter Z: 

dZ 


dt 


= {a + ib)Z - (c id)\Z\^Z + rjz, 


( 5 ) 


where a, b, c(> 0), d are real variables, i = \/^, and rjz is a complex noise term. For a > 0, 
the system starts to oscillate with a mean amplitude Eq. ([S]) can be decomposed 

into two Langevin equations for the amplitude r and the phase 6 oi Z = re*®: 


dr o / \ , 2 / \ 

— = ar - cr^ + r]r{t) , -j-= b - dr +'ne{t), 
dt dt 


( 6 ) 


where rjr and rjo are the white noises of the amplitude and the phase. For simplicity, we 
consider the case where 77 ^ and rjo are uncorrelated with constant strength and Ag re¬ 
spectively. The average phase velocity is u{r) = {dO/dt) = b — dr"^. 

It is clear from Eq. |6]that detailed balance is broken and the system is dissipative. To 
compute the free energy dissipation, we first determine the phase-space probability distri¬ 
bution function P{r,6,t), which follows the Fokker-Planck equation: 

dP 
dt 


* ® V - ^ A] _ k. (T) 


r dr 2 dr do ^ 2 dO ^ r dr dO 

where Jr and Jg are the probability density fluxes in phase space. Since u{r) does not 
depend on 9, the steady state probability distribution Ps{r, 6) only depends on r: 

2(cr"^/4 — ar‘^/2) 


Ps{r, 0) = P{r) = A exp [- 


A^ 


( 8 ) 


where A = [27r f exp [—2(cr'^/4 — ar^/2)/Aj,]rdr] ^ is the normalization constant. From 
Eq. [HI, the flux vanishes in the r-direction Jr = 0. However, there is a finite flux in the 






^-direction Jeir) = uj{r)P{r). We compute the system’s entropy production rate S [371. l38|. 
from which we obtain the minimum free energy dissipation (see SI for details): 


W = kuP 





+ ^ ^^ ]rdrd6 = ksTf, ^ ^ 


A,P AgP 


A, 


e 


(9) 


where Te is an (effective) temperature of the environment, we set ksTe = 1 here. 

The phase diffusion constant is determined by expanding the phase velocity a;(r) around 
r = Vg, the most probable amplitude from P(r). This leads to d6/dt = u{rs) + (36r(t) Prigit), 
with fd = duj{rs)/dr = —2d^^/a/c. The period of the oscillation is T = 27r/uj{rs), and the 
phase fluctuation 66 = 6 — uj{rs)t follows diffusion with the diffusion constant given by: 

De = ^ + A,. (10) 


From Eq. lO MTOl the relation between phase diffusion and energy dissipation emerges: 

where ITc = 0 because = 0, something that is not generally valid (see SI and Fig. S4 
in SI for a more general case of the noisy Stuart-Landau equation). The two constants, 
C = Do = and ITo = (w^)T, depend on the details of the system. 

EqfTTl has the same form as Eq. 0] obtained empirically from studying different biochem¬ 
ical networks. Analysis of the noisy Stuart-Landau equation clearly shows that free energy 
dissipation is used to suppress phase diffusion to increase the coherence of the oscillation. 
Though parameters in this relation may depend on the details of the system, the inverse 
dependence of phase diffusion on energy dissipation appears to be universal. 



III. DISCUSSION 

Oscillations are critical for many biological functions that require accurate time control, 
such as circadian clock, cell cycle, and development. However, biological systems are inher¬ 
ently noisy. The phase of a noisy oscillator fluctuates (diffuses) without bound and eventually 
destroys the coherence (accuracy) of the oscillation. Specihcally, the number of periods W 
in which the oscillation maintains its phase coherence is given by W = Tc/T = aT/D, 
which decreases with the phase diffusion constant. Here, our study shows that free energy 
dissipation can be used to reduce phase diffusion and thus prolong the coherence of the 
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oscillation. A general relationship between the phase diffusion constant and the minimum 
free energy cost, as given in Eq. 01 holds true for all the oscillatory systems we studied here. 
The amplitude fluctuations also decrease with free energy dissipation (see Fig. S5 in SI for 
details), as fluctuations in phase and amplitude are coupled in realistic systems. Our study 
thus establishes a cost-performance tradeoff for noisy biochemical oscillations. 

How do biological systems use their free energy sources (e.g., ATP) to enhance the accu¬ 
racy of the biochemical oscillations? As illustrated in Fig. 5a, a biochemical oscillation can 
be considered as a clock, which goes through a series of time-ordered chemical states (green 
dots) during each period. These chemical states are characterized by the conformational and 
chemical modification (e.g., phosphorylation) states of the key proteins or protein complexes 
in the system. The forward transition from one state to the next is coupled to a PdP cycle 
(blue arrowed circle) driven by hydrolysis of one ATP molecule. For each forward step, 
the reverse transition introduces a large error in the clock. The system suppresses these 
backward transitions by utilizing the ATP hydrolysis free energy. However, this is just one 
half of the story. Even in the absence of the reverse transition, the time duration between 
two consecutive states is highly variable due to the stochastic nature (Poisson process) of 
the chemical transitions. A general strategy of increasing accuracy is averaging [? ]. In 
the case of biochemical oscillations, each period may consist of multiple steps, each powered 
by at least one ATP molecule. As a result of averaging, the error in the period should 
go down as the number of steps increases. Specifically, we expect that the variance of the 
period cr|.(= D/T) should be inversely proportional to the total number of ATP hydrolyzed 
Natp ot T/Pcyc in each period T, where Tcyc is the average PdP cycle time, which is essen¬ 
tially the ATP turnover time. Consequently, the number of coherent period N^. = aT/D 
should be proportional to the number of ATP hydrolyzed in each period. We checked this 
prediction by varying the kinetic rates in the PdP cycle to change Tcyc (see Methods section 
for details). In Fig. 5b, it is shown that the accuracy of the oscillation (clock), as measured 
by Nc, is enhanced by the number of ATP molecules hydrolyzed in each period. This re¬ 
sult reveals a general strategy for oscillatory biochemical networks to enhance their phase 
coherence by coupling to multiple energy consuming cycles in each period. Interestingly, ap¬ 
proximately 15 ATP molecules are consumed per KaiC molecule per period in the circadian 
clock of cyanobacteria (SOj. 

Biological systems need to function robustly against variations in its underlying biochem- 
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ical parameters (rates, concentrations) 1^. For oscillatory networks, the free energy 


dissipation needs to reach a critical value (Wc) to drive the system to oscillate. We showed 
here that additional free energy cost in excess of Wc is needed to make the oscillation more 
accurate, as demonstrated explicitly in Eq. IH In addition to this accuracy-energy tradeoff, 
we found that larger energy dissipation can also enhance the system’s robustness against its 
parameter variations. Take the activator-inhibitor model, for example: the concentrations 
of enzyme E (Et) and phosphatase K (Ex) may vary from cell to cell. We search for the 
existence of oscillation in the {Et, Kt) space for different values of 7 . Robustness is defined 
as the area of the parameter space where oscillation exists. As shown in Fig. S 6 in the SI, 
the robustness increases as the system becomes more irreversible, i.e., when more free energy 
dissipation is dissipated. This suggests a possible general tradeoff between the functional 
robustness and energy dissipation in biological networks. 
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V. METHODS 


Simulation Methods. The Gillespie algorithm ^ is used for the stochastic simulations 
of the reaction kinetics. For given kinetic rates and the volume V, we simulated 1000 
trajectories starting with the same initial condition. For the jth trajectory, we obtained its 
ith peak time tij from the trajectory Xj{t) after smoothing (smooth function in MATLAB 
was used). The peak positions for two trajectories are shown in Fig. 2a. For all the 
trajectories, we computed the mean of their ith peak time rrii = variance 

af = —rrii)^/(N — l), where N is the total number of trajectories. The average period 

T is given by T = rrii/i. Asymptotically, af depends linearly on rrii (Fig. 2d), and the slope 
of this linear dependence is the peak time diffusion constant D, which has the dimension of 
time. The phase diffusion constant is linearly proportional to D\ D^j, = { 271 )“^D/T. For 
the repressilator and the brusselator models, we simulated the stochastic kinetic equations 
to a sufficiently long time (10000 periods) to obtain the time-averaged distribution P{x), 
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where x represents the phase space. We used P{x) to compute free energy dissipation. 

Random Sampling in the {[ATP], [ADP], [PJ) space is performed (in log scale) in the 
region log^o e [2, 5], log^o e [-3,1], log^o G [-3,1] by using Latin hyper¬ 

cube sampling (the Ihsdesign function in MATLAB). The reference concentrations [ATPjo, 
[APPjo, and [Pj]o are set to unity and their actual values are absorbed into the baseline 
reaction rates /-i,o and f- 2 , 0 , which are given in the legend of Fig. 4. 

ATP consumption. In the activator-inhibitor model, the ATP consumption rate is 
Ratp = ^{Jp ~ Jp)i where and J~ are the fluxes for the E ^ Ep and Ep ^ E reactions, 
respectively. We varied the overall reaction kinetics, e.g., Tcyc and the ATP consumption rate, 
by introducing a timescale factor B for all four rates di = d 2 = Bd, fi = f 2 = Bf-, where 
d = 15, / = 15 are the original values used in this paper (see SI). By changing the rates this 
way, the free energy release of ATP hydrolysis AG = — In 7 = ln(ai/ia 2 / 2 /(di/_id 2 /- 2 )) is 
unchanged. We varied B G [0.2,2], and computed the total number of ATP consumed per 
period Natp = Rat pdt and for Fig. 5b. 
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FIG. 1 (preceding page). Different network motifs and the corresponding biochemical oscilla¬ 
tory systems, (a) Illustrations of three network motifs for oscillation: activator-inhibitor, re- 
pressilator, and substrate-depletion, (b) The activator-inhibitor model with a phosphorylation- 
dephosphorylation (PdP) cycle. R and K catalyse two opposing reactions E ^ Ep (phosphoryla¬ 
tion and dephosphorylation) through different intermediate complexes ER and EpK. Ep activates 
both R (activator) and X (inhibitor). X inhibits R by enhancing its degradation. Parameter 
7 = dif-id 2 f- 2 /{aifia 2 f 2 ) is introduced to characterize the reversibility of the system, (c) The 
“repressilator” model of cell cycle in eukaryotic cells. In the simplified network, CDKl activates 
Plkl, Plkl activates APC, and APC degrades CDKl (dashed line), forming the mutually ac- 
tiving/inhibiting loop. Other intermediates are ignored here, (d) The brusselator model with 
detailed reactions. A and B are constant sources, (e) The glycolysis network. The allosteric en¬ 
zyme’s protomer has two states, R (binding with P) and T (unbinding with P), and only R has 
the catalysis activity. Each Rij, with i = 1,2, ••• ,nj and j = 1,2, ••• ,nj represent the num¬ 
ber of S and P bound to R, here we used n* = nj = 2. Each Rij can undergo reactions of 
Rij + S ^i,j + Detailed descriptions and rate values are given in SI. 
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FIG. 2. Correlation and phase diffusion in the activator-inhibitor model with V = 50,7 = 10“^. 
(a) Two noisy oscillation trajectories, with the peaks labeled by circles and squares, (b) Auto¬ 
correlation function (defined in Eq. 1) of the inhibitor X. C{t) decays exponentially with cor¬ 
relation time Tc = 37.7. (c) Raster plot of the peak times for 500 different trajectories starting 
with the same initial condition. The distributions of the peak times for each consecutive peaks are 
shown by red lines. The peak time variance is shown, (d) Peak time variance goes linearly 
with the average peak time, with the linear coefficient defined as the peak time diffusion constant. 
Here, the diffusion constant D = 0.2 and a = TcD/T'^ ~ 0.07. 
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FIG. 3. Relation between the dimensionless diffusion constant (D/T) and free energy dissipation 
per period per volume (AVF, in units of ksT) for the four oscillatory systems. Detailed descriptions 
of the models and parameters can be found in SI. The relationships for different volumes collapse 
onto the same curve when the peak time diffusion constant is scaled by V, as shown in the insets. 
The black dashed line in the activator-inhibitor model indicates the value of AW if we assume 
that hydrolysis of one ATP molecule provides ~ 12kBT energy, which corresponds to 7 ~ 10“^'^. 
All the data can be well fitted with Eq. 4: V x D/T = C + Wq/{W — Wc) (lines in insets), 
where Wc is determined from the critical value 7 c, Wq and C are from fitting. The parameters are: 
(a) activator-inhibitor, Wc = 360.4, VFo = 447.3 ± 55.8, C = 0.28 ± 0.16; (b) repressilator, Wc = 
1.9, Wo = 17.7 ± 3.9, C = 5.5 ± 2.5; (c) brusselator, Wc = 93.1, Wq = 1135 ± 142, C = 0.27 ± 0 . 11 ; 
(d) glycolysis, Wc = 67.4, Wq = 135.8 ± 3.4, C = 0.025 ± 0.019. 
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FIG. 4. The dependence of phase diffusion on the ATP, ADP, and Pi concentrations. We studied 
the activator-inhibitor model with 300 randomly chosen parameters of dimensionless [ATP], [ADP] 
and [Pi] (see Methods). The affected kinetic rates are ai = ai^o[^^.f’])/-i = f-i,o[-^DP], f -2 = 
f- 2 ,o[Pi] with oi^o = 0.1,/_i^o = f- 2,0 = 1- We chose V = 100. (a) D/T versus AW for the 
300 different parameter choices. All the points lie above an envelope curve, which follows Eq. 
4 with D/T = 1.94/(AVF — 400) -|- 0.0036. Points in square indicate the points shown in Fig. 
3. For a given value of D/T, the corresponding minimal energy dissipation AWmin is computed 
according to the fitted envelope curve. The efficiency is defined as E = AWmin/AW. Colors 
of the points indicate the efficiency, (b) Distribution of the 300 randomly sampled points in 
the parameter space {[ATP], [ADP], [Pi]). Colors of the points indicate the efficiency as in (a). 
Points with high efficiency are clustered, (c) Distribution of [ATP], and [ADP]/[Pi] for parameter 
choices with high efficiency E > 0.75. The most probable parameter values for high efficiency are 
[ATP] « 10^, [ADP] « [Pi], which corresponds to ai = 02 ,/-i = f -2 in the kinetic equations. 
This result indicates that high efficiency is achieved when the kinetic rates in the two halves of the 
PdP cycle (phosphorylation and dephosphorylation) are matched. 
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FIG. 5. Experimental evidence from Ref. 


36l | (blue curve) and Ref. 


35l | (Red curve). The two 


experiments measured the oscillation of KaiC phosphorylation in vitro in media with different 
ATP/ADP ratios. The autocorrelation functions were calculated from the original data and fitted 
by a exponential decay cosine function A cos (27rt/T)e“*/'^=, where T is the period, and Tc is the 
correlation time. ln{ATP/ADP) represents the entropic contribution to the free energy, (a) The 
period T is robust against changes in the ATP/ADP ratio, (b) T/tc = a~^D/T decreases with 
\\i{ATP/ADP) and eventually saturates at large ATP/ADP ratio, consistent with our theoretical 
prediction. 
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FIG. 6. Oscillation coherence increases with the number of ATP hydrolyzed per period, (a) 
Illustration of a biochemical oscillation as a clock in phase space. The intermediate states (green 
dots) are represented as the ’’hour ticks” of the clock. The transition from one tick to the next 
is coupled with a ATP hydrolysis cycle. The free energy release AG from the hydrolysis cycle 
powers the forward transition (thick solid arrow) and/or suppresses the backward transition (thin 
dotted arrow). The number of ATP consumed per enzyme molecule in each period T is given by 
T/rcyc, where Tcyc is the average cycle time, (b) The accuracy of the oscillation, characterized by 
the number of correlated (coherent) periods Nc, increases linearly with the total number of ATP 
consumed per period Natp before saturating at very high Natp- We varied Natp hy changing 
the cycle time (see Methods for details), we used V = 100 here. 
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Supplementary Information for 

“The free energy cost of accurate biochemical 
oscillations” by Cao et al. 

I. DESCRIPTIONS OF THE FOUR MODELS 

Here, we describe the mathematical details of the four models of biochemical oscillations 
studied in this paper. For the activator-inhibitor and the glycolysis model, we only give the 
ordinary differential eqnations (ODE’s) to describe the deterministic part of the chemical 
reactions. The actual simulations of the stochastic reactions were done using Gillespie 
algorithm (see Methods in the main text). For the repressilator and brusselator models, we 
use the chemical master equation (or Langevin equation with Poisson noise) and solve the 
corresponding Fokker-Planck equation: 

= -V(FP - DVP) = -VJ, (SI) 

where F is the force vector, and D is the noise matrix. J is the flux vector for each direction. 

In all onr models, the units of the parameters are composed of concentration (c, arbitrary), 
time (f, arbitrary) and volume (U, arbitrary). The molecule number unit cxV represents 
the real counts of a given molecule in the system. For example, if the concentration of 
enzyme E is = 10(c), and the volume is U = lOO(U), then the total number of enzyme 
E in our system is EtV = 1000. 

A. Activator-inhibitor model 

The main components of the model are the activator R and its inhibitor X. R and X are 
linked in a feedback loop through a phosphorylation-dephosphorylation (PdP) cycle (main 
text Fig. lb). R activates the synthesis of both R and X through phosphorylated enzyme 
E, thus forms a positive feedback; at the same time, X degrades R, thus forms a negative 
feedback. The parameter 7 = (iid 2 /-i/- 2 /(aia 2 /i/ 2 ) is introduced to distinguish wether the 
system is in equilibrium (7 = 1) or non-equilibrium (0 < 7 < 1). The kinetics is described 
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by: 


^ = ko[Ep] + hS-k2[X][R] 

ffl = k,[E,] - k,[X] 

® = MEpK] + di[ER] - a,[E]{[R\ - [ER]) - f- 2 [E][K] 

fFm ^^ 2 ) 

^ = a,[E]m - [ER]) + f.,m[R] - [ER]) - (A + d,)[ER] 

^ = h[ER] + d2[EpK] - a2[Ep][K] - f-i[Ep]{[R] - [ER]) 

= a2[Ep][K] + f-2[E][K] - f2[EpK] - d2[EpK] 

with two mass conservation constraints: [E] + [Ep] + [ER] + [EpK] = Et, [EpK] + [K] = Kt, 
where E^ and Kt are the total concentrations of enzyme E and phosphatase K. Each 
term in the equations represents one of the reactions in the main text (see Fig. lb). We 
take symmetric parameters: ko = ki = k^ = l(t~^),k 2 = l{c~^t~^), k 4 = 0.5(t~^),S = 
0A{c),Kt = 1{c),Et = 10(c), ai = 02 = 100(c"A"A,/i = A = A = A = 15(f^),/_i = 
f -2 = y/yaifi/di{c~H~^). The oscillation onset point is at 7 c = 2 x 10“^. Notice that the 
PdP cycle’s reaction rates are much larger than the reactions of synthesis and degradation 
of R and X, so the total R is almost unchanged in the time scale of the PdP cycles. 


B. Repressilator 


We use the simplihed cell cycle model in[l(], where CDK activates Plkl, and Plkl activates 
APC, which inhibits CDK in return. The (deterministic) negative feedback loop kinetics 
are governed by the following ODE’s, with CDK, Plkl, APC concentrations represented by 
x,y,z, respectively: 


dx 

dt 

dy 

dt 

dz 

dt 


ai - jdi- 


yTLl 




0 : 2(1 - y) 

0:3(1 - z ) 


X 




fx dx 


- P2y = fy-dy 


y 


ns 


xr + y^^ 


fd^z = fz- dz 


(S3) 


where A, A are the synthesis and decay rates of each component. We chose ai = 0.1(cf ^), 02 = 
3(f“A)A = 3(cfA>/^2 = l(f = l(^”A)^i = 0.5(c), iC 2 = 0.5(c), iCs = 0.5(c), rii = 
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8,77.2 = 8,773 = 8, and a 3 {t~^) is taken as the control parameter ranges from 1.0 to 3.0. The 
oscillation onset point is = 0.8. 

The full stochastic dynamics is described by the Fokker-Planck equation (Eq. [Ml) with 
the force vector and noise matrix given by: 

F [fx dx, fy dy, fz *^2]) D 2^-^di(ig\^fx T fy T dy, T d^] 


C. Brusselator 


The deterministic equation of brusselator is 


dx 2 

— = a — x + xy, 
dt 


dt 


= b — x^y. 


The Fokker-Planck equation was derived from chemical master equations (CME) 
which gave: 


F 


( a — X + x^y 
b — x^y 



- 1/2 - 2xy + x^/2 
2xy — x^ 12 


(S4) 



(S5a) 


D = 


1 

W 


a + X + x'^y 

9 

—xy 
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—xy 
b + x‘^y 


(S5b) 


where V is the volume of the system. We used b = 0.4, and varied a G [0.12,0.18] in our 
study. Fig. SI gives two examples of probability distribution P{x,t) and fluxes J in the 
brusselator model. 


D. Glycolysis oscillation 

The glycolysis model in our study is taken from , but we have introduced finite reverse 
reaction rates for the catalysis processes to study the free energy dissipation in glycolysis. 
The enzyme PFK are composed of n protomers, and undergoes allosteric regulation by its 
product P. Each protomer exists two states, R, which has catalytic activity for converting 
substrate S to P, and T, which is inactive. Assuming a quasi-equilibrium of the allosteric 
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B. 


states of PFK [3[, the dynamics of S and P can be written as: 


dS 


= Vi 


nD{l + a)"' ^{1 + 6)'^{ka — k'P) 


dt L + (1 + a)^(l + 6*)^ 

dP nD{l + + 9)"'{ka — k'P) 

dt 


(S6) 


-k.P 


L + (1 + aY{i + ey 

where a = {aiS + k'P)/{k + di),0 = a 2 P/d 2 - Parameters were chosen as D = 500(c),n = 
2 , ai = 02 = 10 (c“^t“^), di = d 2 = 10 (f“^),/c = l{t~^),Vi = 0.2{ct~^),ks = 0 . 1 (f“^),L = 
7.5 X 10®. k\c~^t~^) is the reverse reaction rate of P to S. The oscillation onset point 
is /c' = 4 X 10“^. Stochastic simulations were performed for 4 reactions: synthesis of S, 
degradation of P, catalysis of S to P, reverse reaction of P to S. Here, we combined all 
the enzymatic reactions into two reactions since the transitions between different allosteric 
states of the enzyme are much faster than the slow reactions of substrate injection (n,) and 
product removal (kg). 

The dissipation of the system can be directly calculated by summing up the dissipation 
of all the enzymatic reactions: 

^ nD(l + ar-\li-ffnka -k'P) ^ ^ ^ 

L + (1 + a)"(l + 6»)^ k'diP ^ ' 

where Ns quantihes how many molecules of S are catalysed to P. 

II. ENERGY DISSIPATION DETERMINED FROM SOLVING THE 
FOKKER-PLANCK EQUATION 


Consider a general Fokker-Planck equation 

dP{x, t) 


dt 


= -V(FF - DVP) = -VJ 


(S8) 


where F is the force vector, and D is the noise matrix. J is the flux vector for each direction. 
The system’s entropy is 

S(t) = — f P{x,t) In P{x,t)dx (S9) 

n 

The entropy production rate isj^: 


dt 

Integrated by parts: 


dS{t) ^ ^ j "_|_ 1] dx = [[\nP{x,t) + l]VJdT 


dS{t) 

dt 


= — J^VlnP(a;, f)da; 


(SIO) 


(Sll) 
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where is the transposition of J. By definition J = FP — DVP, we have 


V In P = J^D-^F--- 


(S12) 


Finally 



(S13) 


fl) 


The second term is the free energy dissipation rate (also called entropy production rate 
in unit of ksT. 


III. ONSET OF OSCILLATION 


Systems at the onset of oscillation is dissipative. In Fig. S2, we show that at the onset 


of oscillation, i.e., when amplitude is zero, the free energy dissipation is finite positive. The 
free energy dissipation at onset defines Wc that we used in htting main text Fig3. 


IV. EXPERIMENTAL DATA ANALYSIS 


The experimental data were obtained from Ref Gj and Ref[7|. The data were processed 
to calculate the autocorrelation function and htted the autocorrelation function with expo¬ 
nentially decay function A cos(27rf/T)e“^/'^, from which the period T and correlation time r 


could be obtained. See FiefS3l 

V. AMPLITUDE FLUCTUATION AND PHASE DIFFUSION IN THE 

STUART-LANDAU EQUATION 

We first derive the amplitude and phase variance in the simplihed Langevin equations 
from the main text (Eq. 6). The deviation in r can be defined as 


6r = r — r^. 


(S14) 


where = \fajc. Perturbation (5r near follows the equation (in first order approxima¬ 
tion): 



(S15) 
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Following js], as t —oo, we obtain the amplitude fluctuation: 


((5r^(oo)) = 


4a’ 


(S16) 


which only depends on A^. 

The phase is 9{t) = u{r{T))dT. With the expanding expression a;(r(r)) = uj{rs) + 
PSr^r) + Tjoij), we have the phase variance: 


( 02 ) - ( 0)2 = 



0 JO 

t pt 

2 



= /9 

iQ Jo 

Following [ 8 |, as t oo, we have: 


((a;(ri)a;(r2)) - (a;(ri))(a;(r2)))dricir2 

[(5r(ri)5r(r2)) - ((5r(ri))(5r(r2))]dridr2 + Aet. 


(S17) 


(02)-(0)2 = (/32^ + A,)t^D,t 


(S18) 


Clearly, the phase diffusion depends on both Aq and A^. 

Next, we show a general case of the Stuart-Landau equation described by the Langevin 
equations for the two real variables x and y: Z = x + iy. 

^ = (ax - by) - (cx - dy)(x^ + p'O + 

f (S19) 

— = (bx + ay) - (cy + dx){x^ + 0^) + V A 2 r] 2 {t) 

We studied these two equations numerically. We found that if Ai 7 ^ A 2 , then Wc 7 ^ 0. In 
Fig. S4, we show a case with fixed Ai = 0.1 while varying A 2 G [0.05,0.2]. We found that 
the onset, which corresponds to H —)■ 00 , occurs at a finite energy dissipation Wc 7 ^ 0 (Fig. 
S4a). However, the general relationship of D/T = Wq/{AW — Wc) + C, with HA 7 ^ 0 from 
Fig. S4a, holds true (Fig. S4b). 

VI. AMPLITUDE FLUCTUATION 


The amplitude fluctuation can be defined as the dispersion of the stochastic trajectories 
departing from the deterministic trajectory in phase-space: 


02 = / mm{x — Xd)^P{x)dx 


(S20) 


where Xd are the points on the deterministic trajectory. In the four models we studied, the 
amplitude fluctuations decrease with energy dissipation and scale with l/V^, as shown in 
Fig. S5. 
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VII. ROBUSTNESS AND ENERGY DISSIPATION 

For the activator-inhibitor model, the total number of enzyme E and phosphatase K may 
vary in real systems (e.g., from cell to cell). Here we search the parameter space of {Et, Kt) 
in the region Et G [0,1000] and Kt G [0,2], and check wether the system oscillates (with 
amplitude larger than 0.1) for different values of 7. Robustness is dehned as the area in 
the parameter space where oscillation exists. We found robustness increases as the system 
becomes more irreversible or equivalently dissipates more free energy, as shown in Fig. S6. 
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FIG. SI. The dependence of amplitude (^) and period (T) on energy dissipation for the four 
models: (a) activator-inhibitor, (b)repressilator, (c)brusselator, (d)glycolysis. For all the four 
cases, the critical free energy dissipation per period Wc is finite at the onset of the oscillation, i.e., 
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FIG. S2. Probability distributions and fluxes (red arrows) in the brusselator model. (a)a = 0.18, 
near the onset (bifurcation point), (b) a = 0.12, far from the bifurcation point. 
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FIG. S3. Autocorrelation of experimental data. The autocorrelation function were calculated from 
the original data and fitted with Acos(27rt/r)e“*/’^, where T is the period, and r is the correlation 


time. (a). Data from Ref ^ with different ATP ratio. (b).Data from Ref[6| with different ADP 


ratio. 
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FIG. S4. Numerical simulation results of of the general Stuart-Landau equation (Eq. S24) with 
a = 1, c = 1,6 = 2, d = 1, Ai = 0.1. We varied A 2 £ [0.05,0.2]. (a) Relationship between energy 
dissipation AW and noise strength A 2 . When A 2 is large, AW decreases linearly with I/A 2 . The 
dashed line shows when A 2 —> 00 , AW —)• Wc ~ 2967 (black circle), (b) The peak time diffusion 
constant D/T versus AW. The red dashed curve is the htting inverse proportional relation, with 
parameters Wq = 10.2, Wc = 2967,(7 = 0.0011. 
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FIG. S5. Relation between relative amplitude fluctuation bAjA and free energy dissipation AW 
for the four models, (a) activator-inhibitor; (b)repressilator; (c)brusselator; (d)glycolysis. Data for 
different volumes collapses (insets) when we scaled amplitude fluctuation with Ijy/V. 
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FIG. S 6 . The relationship between functional robustness and free energy dissipation, (a) The green, 
blue and red curves in the {Et, Kt) space correspond to the boundaries inside which oscillations 
exist for 7 = 10“^, 7 = 10“^ and 7 = 10“^, respectively. The star indicates the parameters in main 
text Figl a. (b) Robustness, defined as the area of oscillation in the parameter space, increases as 
7 decreases. This means that higher free energy consumptions (on average) is needed for higher 
robustness against parameter variations in achieving oscillatory behaviors. 




